Platelets for Multiscale Analysis in Medical Imaging

This paper describes the development and use of multiscale, platelet-based image reconstruction algorithms in medical imaging. Such algorithms are effective because platelets approximate images in certain (piecewise) smoothness classes significantly more efficiently than sinusoids, wavelets, or wedgelets. Platelet representations are especially well-suited to the analysis of Poisson data, unlike most other multiscale image representations, and they can be rapidly computed. We present a fast, platelet-based maximum penalized likelihood algorithm that encompasses denoising, deblurring, and tomographic reconstruction and its applications to photon-limited imaging. Keywords— Wavelets, Poisson, multiresolution, denoising, tomography I. PHOTON-LIMITED IMAGING Photon-limited imaging arises in several important imaging modalities, including Positron Emission Tomography (PET), Single Photon Emission Computed Tomography (SPECT), and Confocal Microscopy [1]. Frequently diagnostics depend on the quality of these images. Wavelet-based methods are powerful tools for image denoising and have been used in a variety of medical imaging applications; however, most wavelet-based approaches are based on Gaussian approximations to imaging data, and these approximations can be especially inaccurate in the context of photon-limited imaging applications because of the Poissonian nature of the data at low count levels. New multiscale image reconstruction algorithms based on platelets overcome these difficulties. Platelets are localized atoms at various locations, scales and orientations that can produce highly accurate, piecewise linear approximations to images consisting of smooth regions separated by smooth boundaries. Platelet-based techniques for image denoising and reconstruction are demonstrated here in the context of nuclear medicine imaging denoising and confocal microscopy image reconstruction. II. PLATELET APPROXIMATIONS A set of basis functions that approximates or compresses images efficiently is key to reconstructing high-quality images from noisy data. One of the principal ideas underlying wavelets and other multiresolution bases is that of recursive partitions (RPs), which yield sequences of nested image partitions. The wedgelet partition [2] is based on a dyadic, square recursive partition which allows for non-square, “wedgeshaped” partitions only at the final level of the partition. Instead of approximating the image on each cell of the partition by a constant, as is done in a Haar wavelet or wedgelet analysis, we can approximate it with a planar surface. We define R. Nowak was partially supported by the National Science Foundation, grant no. MIP–9701692, the Army Research Office, grant no. DAAD1999-1-0349, the Office of Naval Research, grant no. N00014-00-1-0390. R. Willett was partially supported by the National Science Foundation Graduate Student Fellowship. a platelet fS(x, y) to be a function of the form fS(x, y) = (ASx+BSy + CS) IS(x, y), (1) where AS , BS , CS ∈ IR, S is a dyadic square or wedge associated with a terminal node of an RP, and IS denotes the indicator function on S. Each platelet requires three coefficients, compared with the one coefficient for piecewise constant approximation. In many applications it is beneficial to have the added flexibility of platelet approximations. Image gradients, or smooth transitions between regions of varied intensities, encode information about light emission or reflection as well as surface geometry. In medical imaging applications, radioactive pharmaceuticals diffuse smoothly in regions of homogeneous tissue, resulting in smoothly varying intensities within images of organs. The goal of these applications, however, is often to detect boundaries of or abnormalities within objects under observation. Such regions would be indicated by sharp changes, i.e. edges, in the intensity image. In the context of photon-limited imaging, therefore, it is natural to assume an image model consisting of smooth regions separated by a smooth boundary. The power of platelets is realized in connection with mterm approximations (e.g. approximate images composed of the m best basis functions). As demonstrated in our technical report, for images consisting of smooth regions separated by smooth boundaries, platelet approximations may significantly outperform Fourier, wavelet, or wedgelet approximations [3]. Wavelets and Fourier approximations do not perform well on this class of images due to the boundary. Wavelets adeptly handle point singularities, not curves. Conversely, wedgelets can handle boundaries of this type, but produce piecewise constant approximations and perform poorly in the smoother (but non-constant) regions of images. III. IMAGE DENOISING Platelet-based denoising is possible because we can factor the Poisson likelihood and its parameters into a series of terms based on the image partition. In earlier work, we extended the class of multiscale likelihood factorizations to encompass likelihoods generated by platelet analyses [3]. Multiscale likelihood factorizations are an alternative probabilistic representation (i.e., in addition to that of the original likelihood) of the information in the set of observations x given an underlying intensity array μ [4]. The likelihood factorization allows the likelihood of the entire image to be represented in a tree structure in which both likelihoods and parameter penalties of children are inherited by parents. Using this, it is possible to optimally prune a recursive partition of the data using a fast algorithm reminiscent of dynamic programming and the CART algorithm [3]. The optimal pruning solves a maximum penalized likelihood estimation, wherein the penalization is based on the complexity of the underlying partition. The complexity of a given partition is proportional to the total number of terminal cells. Our goal here is to maximize the penalized likelihood function Lγ(μ) ≡ log p(x | μ) − γ {#θ}, (2) where p(x |μ) denotes a likelihood factorization and {#θ} is the number of parameters in the vector θ (one for each constant terminal cell, three for each platelet terminal cell). The constant γ > 0 is a weight that balances between fidelity to the data (likelihood) and complexity regularization (penalty), which effectively controls the bias-variance trade-off. Solving (P̂, θ̂) ≡ arg max P,θ Lγ(μ(P,θ)) (3) and setting μ̂ ≡ μ(P̂, θ̂) yields a maximum penalized likelihood estimator (MPLE). Larger values of γ produce smoother, less complex estimators; smaller values of γ produce more complicated estimators. We have studied the performance of the MPLE in simulated intensity estimation experiments and investigated the effect of γ; our experiments reveal that γ = 15 log(#counts) consistently results in a low MSE over a broad range of intensity levels. This setting for γ is used in all experiments described in this paper. The algorithm is detailed in our technical report [3], and finds the globally optimal solution to (3), as desired. IV. PHOTON-LIMITED IMAGING APPLICATIONS Platelet-based estimation techniques are here applied to nuclear medicine data for denoising and confocal microscopy data for deblurring. A. Nuclear Medicine Denoising Figure 1 (a) depicts an image of the spine obtained from a nuclear medicine study. Here we applied a technique called “averaging over shifts” or “cycle-spinning”. This entails circularly shifting the raw data by a few pixels, denoising, and then shifting the estimate back to its original position. The result of the application of our denoising algorithm is displayed in Figure 1 (c). Comparisons with the result of hereditary Haar wavelet denoising (in Figure 1 (b)) demonstrate the superiority of platelet representations of image structures. The improvements in edge clarity are most visible in the spinal column and hips. Platelet-based techniques can also be used in tomographic reconstruction applications such as SPECT, as detailed in our technical report [3]. B. Confocal Microscopy Confocal microscopy is used to obtain volume images of small fluorescent objects with high spatial resolution. Due to (a) (b) (c) Fig. 1. Denoising in nuclear medicine. (a) “Raw” nuclear medicine spine image (256× 256 pixels). (b) Haar-based MPLE. (c) Platelet-based MPLE the geometry of these microscopes, a “blurring” is introduced into the measurement process. Since the arrival of fluorescence light at the photo multiplier tube can be modeled as a Poisson process, the “de-blurring” and estimation process may be viewed as a Poisson inverse problem well suited to the application of iterative platelet estimation using the EM algorithm, as detailed in [3]. Figure 2(a) is an image of dendritic branches (with spines) of CA1 pyramidal neurons. The stopped EM-MLE result is displayed in Figure 2(b) and the platelet-based MPLE appears in Figure 2(c). After several iterations the EM-ML estimate worsens considerably with each subsequent iteration. In contrast, the platelet estimate converges, eliminating the need to choose which iteration is the best stopping point, as done with the EM-MLE.